Improved simultaneous multislice magnetic resonance imaging using total variation regularization
Ma Ya-Jun1, Li Sha1, 2, Gao Song1, †
Department of Medical Physics, Health Science Center of Peking University, Beijing 100191, China
Department of Radiotherapy, Beijing Cancer Hospital, Beijing 100142, China

 

† Corresponding author. E-mail: gaoss@pku.edu.cn

Project supported by the National Natural Science Foundation of China (Grant No. 61671026) and the Natural Science Foundation of Beijing, China (Grant No. 7162112).

Abstract

Controlled aliasing in parallel imaging results in higher acceleration (CAIPIRINHA) for simultaneous multislice imaging has been proposed recently, which combines multiband excitation and phase cycling techniques to reduce scan time and improve subsequent imaging reconstruction. In this work, the total variation (TV) regularization method is used to further improve CAIPIRINHA. The TV regularization uses an edge-preserving prior, which establishes a relationship between neighboring pixels for image reconstruction. It reduces artifacts and suppresses noise amplification simultaneously. The results are presented with a standard eight-channel head coil with an acceleration factor of 4, where the TV-regularized CAIPIRINHA generates an improved reconstruction as compared with a typical nonregularized CAIPIRINHA.

1. Introduction

Data acquisition time of magnetic resonance imaging (MRI) is one of the most important considerations for clinical practice. Parallel imaging is a widely used method for scan acceleration; however, it features intrinsic signal-to-noise ratio (SNR) loss due to the reduced data collection.[15] Simultaneous multislice imaging (SMS) is a promising fast imaging technique, which utilizes a multiband radio frequency (RF) pulse to excite multiple slices simultaneously.[68] Because of increased slice-field of view (FOV), the SMS has a better SNR performance as compared with typical parallel imaging techniques. For image reconstruction, the SMS utilizes the inherent spatial sensitivities of multichannel coils to separate the simultaneously excited multiple slices. However, its performance is strongly dependent on the geometry of the coil array. It is difficult to separate slices that have close coil sensitivity profiles, which lead to the SNR and aliasing artifacts decreasing significantly.

The SMS technique has become popular since controlled aliasing in parallel imaging results in higher acceleration (CAIPIRINHA) was proposed in 2005.[9] Using the phase cycling technique,[1012] the simultaneously excited slices shift different spatial distances, thus modifying the appearance of the aliasing artifacts, which reduces their dependence on the coil array geometry. Therefore, the reconstructed images produced by CAIPIRNHA show larger SNR and reduced aliasing artifacts as compared with the images generated from the original parallel imaging method. Recently, a conjugate gradient (CG) iteration reconstruction framework has been developed for the CAIPIRINHA technique.[13] However, the images produced by the CG method without regularization still contain notable aliasing artifacts and noise contamination.

In the recent decades, regularization methods have been widely used in parallel imaging both in image domain and in k-space domain. Two most widely used regularization methods are Tikhonov regularization and total variation (TV) regularization. The Tikhonov regularization has been extensively investigated primarily due to the existence of a closed-form solution. However, it is more suitable for k-space domain.[2] The TV regularization has been widely used in image-domain reconstruction, such as SENSE, because it is an image edge-preserving regularization method for imposing a relationship between neighboring pixels in images.[3,14] At the same time, the noise in the reconstructed images can also be effectively suppressed.

In this study, the TV regularization method is introduced and applied to CAIPIRNHA reconstruction. In-vivo data results demonstrate that TV regularization is a powerful technique for suppressing noise, maintaining fine imaging details, and reducing aliasing artifacts.

2. Theory
2.1. Conventional multislice CAIPIRNHA

The CAIPIRNHA technique modifies the appearance of the simultaneously excited multislice images by phase cycling technique. It can be easily implemented using designed multiband RF pulses with different RF phase modulations between all the excited slices. Thus, each excited slice shifts different distances in the image domain.

Without loss of generality, we assume a two-fold scan time reduction (i.e., simultaneously excite two slices) with an N-coil receiver array. For the first slice, the RF phases are kept identical and its j channel image is the discrete inverse Fourier transformation of the underlying discrete k-space signal S1, j as where Δkx = 2π/FOVxky = 2π/FOVy) represents the distance between M(N) equidistantly sampled k-space points in the x(y) encoding direction, and index m (n) is the actual encoding step.

For the second slice, the RF phases increase linearly. Thus, like the well-known Fourier shift theorem, a spatial shift of Δy in the phase encoding direction can be expressed as

Then, the folded image acquired in practice can be expressed as the summation of Eqs. (1) and (2) as where C1, j and C2, j are the coil sensitivity profiles of the two slices, and ρ1 and ρ2 are the spin densities of the two slices, respectively. The spatial shift Δy changed the folded positions of the simultaneously excited two slices. It can increase difference between the two coil sensitivities at the same pixels within the folded image. This modification is beneficial to the unfolded procedure, such as SENSE. If we know the coil sensitivity profile, the two-slice image can be reconstructed using the CG iteration reconstruction algorithm. In addition, the SNR of the final image can be expressed as where NS is the number of slices, g is the geometry factor, and R is the acceleration factor. The resulting SNR of the CAIPIRNHA is at least xs times that of the traditional parallel imaging with the same acceleration factor.

2.2. TV-regularized MS CAIPIRNHA

The images reconstructed from the traditional unfolded algorithm (CG) still contain notable noise and aliasing artifacts when large acceleration factor is needed. TV regularization is introduced for better eliminating the aliasing artifacts and suppressing the noise amplification.

A TV-regularized reconstruction framework is proposed as follows: where f is the unfolded image and d is the acquired k-space data from all the channels, E is the spatial encoding matrix combined with gradient and coil sensitivity encodings and the matrix element is as follows: where Cns,j is the coil sensitivity profile of the NS slice and j channel. This presentation is based on forward Fourier transformation, while in the previous section the inverse Fourier transformation was used. However, they are equivalent to each other with the same physical nature.

The second term on the right-hand in Eq. (5) is the TV regularization term, which is defined as a function of the image gradient where ∇x and ∇y denote the gradients along the x and y directions in the 2D image, respectively, and | | denotes the complex modulus. λ is the TV regularization parameter.

Using the TV regularization term is based on the fact that highly oscillating noise and irregular folded artifacts often increase the TV value of an image. The λ is the parameter to find the balance between the first least-square norm and the second TV term, and its value can be estimated by the L-curve method,[15,16] or can be simply generated from the automatic search method.

3. Materials and methods

For identifying the effects of TV regularization on the CAIPIRNHA, in-vivo data were acquired from human brain of a healthy volunteer on a Philips 3T MRI scanner using a normal spin-echo sequence. Written informed consent approved by the Institutional Review Board was obtained prior to his participation. Imaging parameters were as follows: TR/TE = 490 ms/11.6 ms, matrix size = 256 × 196, field of view = 230 mm × 198 mm, number of slices = 18, slice thickness = 4 mm, and slice gap = 1 mm. The imaging data were acquired with an eight-channel head coil.

Simultaneously excited data were derived from the summation of four slices (these slices were chosen from the acquired 18 slices: the 1st, 3rd, 8th, and 10th slices) with different phase cycles: 0, π/2, π, 3π/2. Figure 1 shows two channel images.

Fig. 1. (a) and (b) Two channel images of the simultaneously excited image-domain data with an acceleration factor R = 4.

Coil sensitivity profiles of the four excited slices were obtained by the self-calibration method using 48 phase-encoding lines at the center of their k-space.[15] Then, a fixed-point iteration algorithm was used to solve the nonlinear minimization problem of Eq. (5). All the data processing was implemented in MATLAB 2012a (64 bit; Mathworks, Natick, Massachusetts, USA) on an i7 type Intel CPU. In addition, CG iteration reconstruction was also implemented for comparison.

For quantitative comparison, the normalized mean squared error (NMSE) is defined as the normalized square of the difference between the reconstructed image (Irec) and the reference standard image (ISOS: it was obtained by the sum-of-square (SOS) method for each single acquired slice). Therefore, the NMSE is expressed as follows:

4. Results and discussion

The first three rows in Fig. 2 show four slices produced by the SOS, CG, and TV reconstruction methods. The same regions in all the images are amplified for better comparison. Noise can be seen in the images by CG reconstruction, and the images using TV-based reconstruction have almost the same image quality as that of the standard SOS images.

Fig. 2. (color online) Imaging results produced by SOS (1st row), CG (2nd row), and TV (3rd row) reconstruction methods, and error maps of CG (4th row) and TV (5th row) reconstructed images and their values are normalized and displayed with the same color scale.

The error maps calculated by the two equations: ErrCG = |ICGISOS| and ErrTV = |ITVISOS| are shown in the fourth and fifth rows in Fig. 2. They are normalized and displayed with the same color scale. The intensities of noise and artifacts are smaller and lower in the TV images than in the CG images.

The NMSE values are shown in Fig. 3, showing the total errors of the TV and CG images. It can be seen that the NMSE values of the TV images are always lower than those of the CG images.

Fig. 3. (color online) NMSE lines of the CG and TV methods.

In this study, the simulated multislice data are used for image reconstruction. Though the simulated data may be slightly different from the real data, it still provides important information for comparing the two reconstruction methods. Besides, it is easy to quantitatively compare the two reconstruction methods since the standard images are known from the simulated data. In our further studies, we will compare the two reconstruction methods for the real data.

The CAIPIRNHA technique has largely improved the image quality compared with the traditional parallel imaging.[2,9] In this work, we find that the image quality of TV-regularized reconstruction for CAIPIRNHA is much better than that of CG reconstruction. However, equation (5) is a nonlinear minimization problem that requires a longer computation time to be solved by the TV reconstruction than by the CG method. The CG reconstruction needs 14.7 s for 80 iterations and the TV reconstruction needs 54.5 s for 20 iterations. To realize real-time calculation in our future work, advanced tools, such as Compute Unified Device Architecture (CUDA) programing based on Graphics Processing Unit (GPU), can be used.

Initial values are crucial for the non-linear TV reconstruction. In this study, we employ the zeros-filled data with 48 center k-space lines, which is used for coil sensitivity calculation, as the initial images for reconstruction. With this initial image, less iterations are needed and more accurate results are obtained.

In addition, the selection of a proper TV regularization parameter λ still needs further investigating.[11] Here, we employ an automatic search method to find the best TV regularization parameter that is accurate but needs time.

To further accelerate the MR data acquisition, the CAIPIRNHA technique can be combined with the typical parallel imaging method. The proposed TV regularization can also be applied to this case.

In addition, more data with different image contrasts can be employed to further test the TV-regularized CAIPIRNHA.

5. Conclusions

In-vivo data results demonstrate that our proposed TV regularization method is significant in the reconstruction of simultaneous multislice images with improved SNR and reduced aliasing artifacts as compared with the traditional nonregularized method.

Reference
[1] Larkman D J Hajnal J V Herlihy A H Coutts G A Young I R Ehnholm G 2001 J. Magn. Reson. Imaging 13 313
[2] Griswold M A Jakob P M Heidemann R M Nittka M Jellus V Wang J Kiefer B Haase A 2002 Magn. Reson. Med. 47 1202
[3] Pruessmann K P Weiger M Scheidegger M B Boesiger P 1999 Magn. Reson. Med. 42 952
[4] Li S Wang L Zhu Y C Yang J Xie Y Q Fu N Wang Y Gao S 2016 Chin. Phys. 25 128703
[5] Bao S L Du J Gao S 2013 Acta Phys. Sin. 62 088701 (in Chinese)
[6] Ma Y J Liu W Tang X Gao J H 2014 Magn. Reson. Med. 74 217
[7] Fang S Wu W C Ying K Guo H 2013 Acta Phys. Sin. 62 48702 (in Chinese)
[8] Wang L Y Liu H K Li L Yan B Zhang H M Cai A L Chen J L Hu G E 2014 Acta Phys. Sin. 63 208702 (in Chinese)
[9] Breuer F A Blaimer M Heidemann R M Mueller M F Griswold M A Jakob P M 2005 Magn. Reson. Med. 53 684
[10] Glover G H 1991 J. Magn. Reson. Imaging. 1 457
[11] Zu Z L Zhou K Zhang S G Gao S Bao S L 2008 Chin. Phys. 17 328
[12] Zhu Y C Spincemaille P Liu J Li S Nguyen T D Prince M R Xie Y Q Wang Y 2017 Chin. Phys. 26 018701
[13] Yutzy S R Seiberlich N Duerk J L Griswold M A 2011 Magn. Reson. Med. 65 1630
[14] Block K T Uecker M Frahm J 2007 Magn. Reson. Med. 57 1086
[15] McKenzie C A Yeh E N Ohliger M A Price M D Sodickson D K 2002 Magn. Reson. Med. 47 529
[16] Li S P Wang L Y Yan B Li L Liu Y J 2012 Chin. Phys. 21 108703